Transcriptomic analysis of Anopheles gambiae from Benin reveals overexpression of salivary and cuticular proteins associated with cross-resistance to pyrethroids and organophosphates

Background Insecticide resistance (IR) is one of the major threats to malaria vector control programs in endemic countries. However, the mechanisms underlying IR are poorly understood. Thus, investigating gene expression patterns related to IR can offer important insights into the molecular basis of IR in mosquitoes. In this study, RNA-Seq was used to characterize gene expression in Anopheles gambiae surviving exposure to pyrethroids (deltamethrin, alphacypermethrin) and an organophosphate (pirimiphos-methyl). Results Larvae of An. gambiae s.s. collected from Bassila and Djougou in Benin were reared to adulthood and phenotyped for IR using a modified CDC intensity bottle bioassay. The results showed that mosquitoes from Djougou were more resistant to pyrethroids (5X deltamethrin: 51.7% mortality; 2X alphacypermethrin: 47.4%) than Bassila (1X deltamethrin: 70.7%; 1X alphacypermethrin: 77.7%), while the latter were more resistant to pirimiphos-methyl (1.5X: 48.3% in Bassila and 1X: 21.5% in Djougou). RNA-seq was then conducted on resistant mosquitoes, non-exposed mosquitoes from the same locations and the laboratory-susceptible An. gambiae s.s. Kisumu strain. The results showed overexpression of detoxification genes, including cytochrome P450s (CYP12F2, CYP12F3, CYP4H15, CYP4H17, CYP6Z3, CYP9K1, CYP4G16, and CYP4D17), carboxylesterase genes (COEJHE5E, COE22933) and glutathione S-transferases (GSTE2 and GSTMS3) in all three resistant mosquito groups analyzed. Genes encoding cuticular proteins (CPR130, CPR10, CPR15, CPR16, CPR127, CPAP3-C, CPAP3-B, and CPR76) were also overexpressed in all the resistant groups, indicating their potential role in cross resistance in An. gambiae. Salivary gland protein genes related to ‘salivary cysteine-rich peptide’ and ‘salivary secreted mucin 3’ were also over-expressed and shared across all resistant groups. Conclusion Our results suggest that in addition to metabolic enzymes, cuticular and salivary gland proteins could play an important role in cross-resistance to multiple classes of insecticides in Benin. These genes warrant further investigation to validate their functional role in An. gambiae resistance to insecticides. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10261-x.


Background
Insecticide-based vector control approaches such as insecticide treated nets (ITNs) and indoor residual spraying (IRS) are core methods to break humanvector contact, thus reducing transmission of malaria [1].The insecticides traditionally used in vector control fall within 4 main classes: pyrethroids, carbamates, organochlorines, and organophosphates, with pyrroles and neonicotinoids being recently repurposed from agricultural pesticides [2].Because of this reliance on insecticide-based control methods, the emergence and spread of insecticide resistance in mosquito vectors is becoming a serious threat.
Insecticide resistance is caused by multiple mechanisms including target site insensitivity, such as knockdown resistance (kdr), and metabolic resistance brought about by the increased production of enzymes capable of breaking down insecticides [3].Knockdown resistance (kdr) arises from mutations in the voltage-gated sodium channel gene which is a target of pyrethroid and organochlorine insecticides, while the acetylcholinesterase (ace-1) gene is a target of organophosphate and carbamate insecticides.Mutations such as L1014S and L1014F in the kdr gene [4] and the G280S mutation in the ace-1 gene [5] cause structural changes that decrease the ability of certain insecticides to bind with their target sites.Additionally, the overexpression of genes encoding detoxi cation enzymes (e.g., cytochrome P450s, glutathione S-transferases and carboxylesterases) are known to contribute to metabolic resistance in mosquitoes [6,7] and can be markers of IR in mosquito populations [8].Because these changes resulting in IR have a genetic basis, robust sets of molecular markers could provide sensitive and timely diagnoses of insecticide resistance in mosquito populations, thereby enabling the implementation of effective insecticide resistance management (IRM) strategies.
Next-generation sequencing approaches such as RNA-seq and whole genome sequencing (WGS) enable the understanding of the genetic variations that result in changes to mosquito biology and behavior in response to environmental factors, including insecticide pressure.Transcriptomic analyses allow for the investigation of gene expression and polymorphic variations associated with speci c phenotypes [9,10].Furthermore, it allows the identi cation of novel candidate genes related to a speci c insecticide or multi-insecticide resistance [9,11,12].A recent study of the malaria vector Anopheles arabiensis found multiple highly overexpressed genes related to cuticular-associated proteins and salivary gland proteins associated with pyrethroid and organophosphate resistance, suggesting roles of these lesser-understood gene groups in cross-resistance [12].
The objective of this study was to use transcriptomic data to identify candidate genes associated with insecticide resistance in An. gambiae s.s collected from Djougou and Bassila, two sites in Benin with different levels of resistance to pyrethroid or organophosphate insecticides.

Results
Resistance pro le of An. gambiae from Bassila and Djougou The mortality of the mosquitoes in each insecticide bioassay is presented in Fig. 1 and Additional le 1.There was a signi cant difference in the mortality rate between mosquitoes from Bassila and mosquitoes from Djougou (P < 0.05).Mosquitoes from Djougou were more resistant to alphacypermethrin and deltamethrin than mosquitoes from Bassila.At 1X exposure to alphacypermethrin and deltamethrin, the average percent mortality in Djougou was zero, while it was 77.5% (sd: 2.28) for alphacypermethrin and 65.5% (sd: 3.03) for deltamethrin in Bassila.At 2X exposure, the average percent mortality in Djougou was 46.7% (sd: 1.06) and 26.68% (sd: 2.5) for alphacypermethrin and deltamethrin, respectively, and in Bassila, it was 82.8% (sd: 0.1) and 80% (sd: 5.54), respectively.At 5X exposure, the average percent mortality in Djougou was 77.4% (sd: 1.65) and 50.8% (sd: 2.98) for alphacypermethrin and deltamethrin, respectively, and 100% and 86.8% (sd: 0.72) in Bassila (Additional le 1).

RNA-Seq data quality control and mapping
Raw reads generated ranged from 45-112 million for mosquitoes from Bassila, 51-108 million for mosquitoes from Djougou and 77-102 million for the susceptible Kisumu strain mosquitoes.After ltering, more than 98% of the reads were retained in all experiments and mapped to the whole reference genome of Anopheles gambiae PEST (VectorBase release 48).The percentage of reads mapped to the reference genome ranged between 61% and 70.6% for Bassila, 51.4% and 76% for Djougou and 69.2% and 75% for the susceptible Kisumu strain (Additional le 2), and 70% and 79% of the alignments (read pairs) were successfully assigned to the exonic features of the gene set AgamP4 (Additional le 3).

Differential gene expression analysis
EdgeR was used to perform differential gene expression (DGE) analysis between the different mosquito groups (R-S), (C-S) and (R-C).A fold-change (FC) > 2 and a false discovery rate (FDR) < 0.01 were used to identify differentially expressed genes (Fig. 2).The DGE results are summarized in Table 1 and Additional le 4. Genes associated with alphacypermethrin resistance were derived from three different comparisons using mosquitoes from Djougou.A total of 1274 genes (698 upregulated and 576 downregulated), 1605 genes (850 upregulated and 755 downregulated) and 25 genes (13 upregulated and 12 downregulated) were signi cantly differentially expressed in the C-S, R-S and R-C comparisons, respectively (Table 1).A total of 4 key upregulated genes were shared between the three comparison sets, with only one characterized protein belonging to the cuticular protein RR-1 family 75 (Additional le 5-A).Comparing R-S and C-S comparisons, genes overexpressed in both sets are shown in (Additional le 6).Interestingly, a carboxylesterase (COE22933), some cuticular proteins (CPR76, CPR75, CPR16, TWDL1, CPR81, and CPAP3-A1b), some members of the cytochrome P450 family (CYP4D17, and CYP12F2) and some salivary gland proteins (salivary secreted peptide (AGAP013060), kDa salivary (AGAP004316)) showed higher fold changes in the R-S comparison compared to the C-S comparison.Along with these genes, DE-cadherin-like isoform X1 (AGAP029696), COX2, TEP1, ESP, methyltransferase 22 isoform X1 (AGAP011879), micro bril-associated glyco 4-like (AGAP010531) and BTB/POZ domain-containing protein KCTD16 (AGAP000805) also showed higher fold changes in the R-S set of genes (Additional le 6).Importantly, a comparison between the R-S and R-C groups available in (Additional le 7) showed that the cuticular protein (CPR75) and the carboxylesterases (COEunkn) were differentially expressed and shared among both comparisons.
Comparing the R-S and R-C groups, CYP6Z3 in Bassila, and CPR9, CPR144 and CYP4H24 in Djougou were shared among both comparisons.In addition, the angiopoietin-like salivary protein (AGAP007041) in Bassila and CYP314A1 in Djougou were overexpressed in only the R-C group, together with other genes such as brinogen A (AGAP011228), CLIPB12, PPO6, and PGRPLB in Djougou (Additional le 7).
Comparing the R-S and R-C groups, no metabolic genes were shared among the two comparisons in Bassila, while CYP6Z3 and GSTD11 were shared in Djougou.Additionally, two carboxylesterases (AGAP011507, COEAE2E), PPO6, PPO9, PGRPLB, and brinogen A (AGAP011228) were overexpressed in only the R-C group in Djougou (Additional le 7).

Genes associated with multiple insecticide resistance
A total of 500 DEGs (265 upregulated and 235 downregulated) were shared by mosquitoes that survived exposure to either alphacypermethrin, deltamethrin or pirimiphos-methyl (Fig. 5-A, Additional le 9).Among those genes, 10 were cuticular proteins (CPAP3-B, CPAP3-E, CPR10, CPR130, CPR15, CPR16, CPR30, CPR75, CPR76, and CPR81), 8 were cytochrome P450s (CYP12F2, CYP307A1, CYP4C26, CYP4C27, CYP6Z3, CYP9J3, CYP9K1, and CYP9L3), 3 were glutathione S-transferases (GSTD7, GSTE2, and GSTMS3) and 3 were unde ned salivary gland protein genes Fig. 5-C, Table 2).Non-synonymous target site mutations Target site mutations in the Ace-1 and in the kdr genes were identi ed through the analysis of single nucleotide polymorphisms (SNPs).The analyzed samples were composed of pools of seven mosquitoes, and a single population included three replicates (pools).The frequency at which polymorphisms appeared in a population was derived from the depth coverage at each position with the contribution of each mosquito in the pools.We found the G280S mutant allele frequency in the Ace-1 gene to be 50% for Bassila compared to 66% for Djougou (Fig. 6).For the kdr gene, we observed the L995F mutation in 33% of samples from Bassila and 100% of samples from Djougou.Other mutations observed in that gene are shown in Fig. 6.

Gene ontology annotation and enrichment analysis
Gene ontology enrichment (GOE) analysis was conducted on differentially expressed genes (up and downregulated) for all R-S comparisons.Gene ontologies are classi ed into three classes: biological process (BP), molecular function (MF) and cellular component (CC) (Additional le 10).
Among BP GO terms, terms related to proton transport (GO:0015986), aerobic respiration (GO:0009060) and mitochondrial electron transport (GO:0006123, GO:0006120) were shared in the BD vs KIS, BP vs KIS, DD vs KIS and DP vs KIS comparisons.Furthermore, GO terms associated with glycolytic processes (GO:0006096) were enriched in BP vs KIS and DP vs KIS; while the carbohydrate metabolic processes (GO:0005975) were enriched in BD vs KIS and DD vs KIS (Additional le 11).Considering the CC GO terms, terms associated with mitochondrial respiratory chain complex I (GO:0005747) III (GO:0005750) and IV(GO:0005751) were enriched in BD vs KIS, BP vs KIS, DD vs KIS and DP vs KIS with mitochondrial proton-transporting ATP synthase complex, and coupling factor F(o) (GO:0000276) was enriched in only the BP vs KIS and DP vs KIS comparisons (Additional le 11).Concerning terms related to MF, macromolecular complex binding related terms (GO:0044877) and sulfur cluster binding (GO:0051539) were enriched in all R-S comparisons except DA vs KIS.Hydrogen ion transmembrane transporter activity (GO:0015078) was enriched in both BP and DP vs KIS, highlighting its strong correlation with exposure to pirimiphos-methyl (Additional le 11).

RNA-Seq data validation using quantitative PCR
The expression patterns of four genes (SG7, CYP9K1, CYP6P3 and COEJHE5E) were validated in relation to two housekeeping genes (40S ribosomal protein S7; RPS7 and Actin5c) (Fig. 7).Most of the qPCR results supported the directionality of the expression level changes observed after RNA sequencing.
Nevertheless, for the COEJHE5E gene, while downexpression was observed in DP vs KIS after qPCR, an overexpression was observed in the RNA-seq dataset.

Discussion
In this study, we analyzed the gene expression pro les of two An.gambiae populations following exposure to key insecticides used for malaria vector control.In addition to detecting differential expression of genes encoding detoxi cation enzymes, genes coding for salivary gland and cuticular proteins were also overexpressed in mosquitoes from both study sites.
The two mosquito populations displayed variable levels of resistance to alphacypermethrin, deltamethrin and pirimiphos-methyl.A higher level of resistance to both alphacypermethrin and deltamethrin was observed in Djougou than in Bassila despite both being located in Donga district in northern Benin.Such heterogeneity in the levels of insecticide resistance despite relative geographic proximity has also been described elsewhere [13,14].The high pyrethroid resistance detected in the Djougou population could be related to the intensive use of pyrethroids in vector control by the NMCP which has intensi ed ITN distribution every three years since 2008 in Djougou [15].In addition, there is intensive use of carbamates, organophosphates and pyrethroid pesticides for agricultural activity in Donga district [15,16] .In contrast, in Bassila, apart from being an agricultural site, no NMCP interventions have been distributed to date.The lower alphacypermethrin and deltamethrin resistance in mosquitoes from this locality might thus be related to reduced selective pressure.For this reason, the increased resistance to pirimiphos-methyl is surprising as it was higher than that of Djougou.
Cytochrome P450s are usually present in abundance in insects [17] and their primary function is to metabolize pheromones and xenobiotics [17][18][19][20].In the Djougou mosquitoes that survived alphacypermethrin exposure, two members of this family, CYP4H17 and CYP4D17, were overexpressed.Both genes have been found to be overexpressed in mosquitoes resistant to permethrin [21,22] with an overexpression of CYP4D17 in mosquitoes exposed to 0.05% deltamethrin [23].Their over-expression only in mosquitoes resistant to alphacypermethrin here is quite interesting given that they were not overexpressed in deltamethrin resistant from the same population.Further, CYP4G16 and CYP12F3 were overexpressed in Djougou mosquitoes that survived both alphacypermethrin and deltamethrin exposure.Overexpression of CYP4G16 has been demonstrated to be important in cuticular hydrocarbon content enrichment and thus reducing the uptake of pyrethroid insecticides [23,24].On the other hand, over-expression of CYP12F3 has not yet been associated with insecticide resistance.Functional validation of these four genes associated with the high intensity of resistance recorded to these class two pyrethroids would con rm the role played in conferring resistance.
Several CYP450s were found to be commonly shared among all R-S comparisons including CYP9K1, CYP6Z3 and CYP12F2.All three were shown to be consistently overexpressed in both pyrethroid and pirimiphos-methyl resistant mosquitoes.CYP9K1 has previously been shown to be involved in the metabolism of deltamethrin and pyriproxyfen in An. gambiae [23,25] and therefore its overexpression in the pyrethroid-resistant groups is not surprising .
Additionally, a similar study in An. arabiensis showed its association with resistance to both organophosphates and pyrethroids [12], consistent with our ndings here of its association with resistance to both chemical classes.CYP6Z3 has previously been associated with pyrethroid metabolism [24] but here we observed its overexpression in survivors exposed to both pyrethroids and pirimiphos-methyl indicating a potential role in cross-resistance as has been reported before [26][27][28][29].CYP12F2 has been previously associated with permethrin resistance in An. arabiensis [30] [31] and was associated with both pyrethroid and organophosphate resistance in our study, suggesting a role in the detoxi cation of both insecticide classes.Moreover, CYP6Z2, CYP9M1 and CYP4H15 were found to be overexpressed in 4 out of 5 R-S comparison groups, with the exception of BD vs KIS.The Bassila population was not as intensely resistant to pyrethroids as Djougou, which suggests that these three genes could potentially contribute to intensi ed pyrethroid resistance as well as pirimiphos-methyl resistance.
Other metabolic genes including three glutathione-S-transferase genes (GSTD7, GSTE2, GSTE4) and two carboxylesterase genes (COEJHE5E, COE22933) were overexpressed in all R-S comparisons.This is consistent with the understanding that GSTs contribute to resistance against multiple insecticide classes [32] .They participate in the detoxi cation of xenobiotics and metabolize secondary products from other metabolic activities (CYP450s, COEs) [32,33] .COEJHE5E has been recently found to be overexpressed in deltamethrin and permethrin resistant An. coluzzii population [23], while genomic signals of resistance to deltamethrin were found around the COE22933 locus in An. gambiae and An.coluzzii [34].A functional validation of these two genes might enhance understanding of their role in pyrethroid resistance.
In addition to metabolic genes, salivary gland and cuticular proteins were also differentially expressed.In our study, a wide variety of genes encoding cuticular proteins (CPs) were found to be overexpressed in R-S comparisons.Among them, CPR10, CPAP3-E, CPAP3-B, CPRd30, CPR130, CPR15, CPR16, CPR76, CPR75, and CPR81 were overexpressed in all R-S comparisons.Their mechanism of function is not well understood but alterations to the mosquito cuticle by thickening or hardening could inhibit the penetration of insecticides and other toxicants [35].Additionally, recent studies highlighted the massive production of chitin to enhance metabolic gene functions, such as CYPs, which could lead to cuticle hardening in Aedes aegypti [36].If this is the case, this mechanism would lead to cross-resistance of insecticides regardless of class.CPAP3s belong to a group of proteins that maintain the structural integrity of the cuticle [37] .CPAP3s with an obstructor-E function have been found to play a role in multiple insecticide resistance in An. gambiae [38].In this study, CPAP3-B was overexpressed in all R-S comparisons making this gene a potential candidate marker for cross-resistance.Other cuticular proteins that belong to the CPR family have been associated with multiple-resistance phenotypes [35] .In this study, CPR10, CPR30, CPR130, CPR15, CPR16, CPR75, CPR81 and CPR76 were overexpressed, although their roles in insecticide resistance have not yet been investigated.Additional cuticular genes were observed to be more speci cally overexpressed in response to exposure to certain classes of insecticides; CPR9, CPR59, CPLCP3, TWDL1, CPAP3-A1b and CPAP3-A1c in response to pyrethroids and CPR144 in response to pirimiphos-methyl exposure.The involvement of CPLCP3 in cuticle barrier formation has been previously described and its overexpression has been associated with mosquito survival following deltamethrin exposure [39].The overexpression of CPAP3-A1b and CPAP3-A1c has also been associated with resistance to permethrin, deltamethrin or lambdacyhalothrin [40][41][42].
Recently, salivary gland proteins have been reported to play a potential role in insecticide resistance [12].In this study, three salivary gland genes stood out among the overexpressed genes in resistant mosquitoes.A threonine serine-rich mucin (SG9: AGAP013423, FC ranging from 2.07 to 2.62) and a salivary cysteine-rich peptide (AGAP011460, FC ranging from 2.03 to 3.38) were overexpressed in mosquitoes resistant to pirimiphos-methyl while salivary secreted mucin 3 (AGAP009473, FC ranging from 2.04 to 2.21) was overexpressed in Djougou mosquitoes resistant to alphacypermethrin and deltamethrin.The salivary protein SG9 is a mucin protein.The role of mucin protein in mosquitoes remains uncharacterized to date, unlike in vertebrates (including humans).They are found in the peritrophic membrane and have been reported to only be induced after blood ingestion in female mosquitoes [43] .Here, we found a slight overexpression (FC ranging from 2.51 to 2.62) of this gene in mosquitoes resistant to pirimiphos-methyl.It has been demonstrated that salivary gland mucin-like protein in Drosophila could perform an immune defense reaction [44], and for this reason, we hypothesize that SG9 could perform a similar function in protecting and defending the respiratory wall against the penetration of organophosphate molecules.In humans, for example, exposure to lowlevel OPs triggers increased mucin secretion in asthmatic patients [45] .This might also be the case in mosquitoes, highlighting the need to further investigate mucin-like protein functions.
On the other hand, AGAP011460 and AGAP009473 are two uncharacterized genes that were both highly overexpressed in mosquitoes resistant to pirimiphos-methyl and both pyrethroids.Further characterization and functional annotation are needed to understand their involvement in insecticide resistance.
In addition to differential gene expression, target site mutations are also important genetic markers of phenotypic insecticide resistance.Mutations such as L995F or L995S in the kdr gene [4] and G280S mutations in the Ace-1 gene [5] have been reported to confer resistance to insecticides.The data on target site mutations presented here are relatively low resolution, as they were extrapolated from RNA-Seq data on pooled mosquitoes.Nevertheless, these data provide some insight into the presence of these mutations in the two populations.The mutation L995F was found in Djougou with an allele frequency of one (100%) and in Bassila with an allele frequency of 0.33 (33%).This is likely the result of evolution due to the selective pressure caused by the intensive use of pyrethroid insecticides in Djougou [15,16] .A frequency of 0.66 (66%) of the G280S mutation which can confer resistance to organophosphates was detected in Djougou.although this population was less resistant to pirimiphos-methyl than Bassila where the mutation frequency was 0.5 (50%).Previous research has shown that the G280S mutation may not be a strong predictor of resistance to all organophosphates [46], which supports this study's ndings.
The overexpression of other genes described above could be important contributors to the resistant phenotypes observed.

Conclusions
The analysis of RNA-Seq data allowed us to describe differentially expressed genes and target site mutations that are associated with pyrethroid and organophosphate resistance in An. gambiae from two locations in Benin: Bassila and Djougou.Multiple genes, including members of the cytochrome P450 family (CYP4H17, CYP4D17 and CYP12F2), salivary gland proteins (SG9, AGAP011460 and AGAP009473) and cuticular proteins (CPR30, CPR130, CPR15, CPR16, CPR76, CPAP3-A1b and CPAP3-A1c) were found overexpressed in resistant mosquitoes after exposure to alphacypermethrin, deltamethrin, and pirimiphos-methyl.The DEGs described here are potential molecular markers of insecticide resistance that could be incorporated into Benin's NMCP insecticide resistance surveillance and management strategy once validated.
Anopheles gambiae from the insecticide susceptible reference Kisumu laboratory strain were reared in the insectary at the Centers for Disease Control and Prevention (CDC), Atlanta, Georgia, USA.Mosquitoes were maintained at a constant 27 ± 2°C and 70 ± 10% humidity on a 14 h:10 h hour light:dark cycle and adults were provided 10% sucrose ad libitum.

Insecticide resistance intensity bioassays
Insecticide resistance intensity assays using modi ed CDC bottle bioassays were conducted on the reared adult mosquitoes using 1X, 2X, and 5X the diagnostic doses of deltamethrin and alpha-cypermethrin and 0.5X, 1X and 1.5X the diagnostic doses of pirimiphos-methyl.Stock solutions (10X) of the diagnostic doses (1X for deltamethrin: 12.5 µg/bottle, alphacypermethrin: 12.5 µg/bottle and pirimiphos-methyl: 20 µg/bottle) were prepared by diluting technical grade insecticide in 50 mL of absolute ethanol.Lower doses were obtained by serial dilution.
For each concentration, four bottles were coated with 1 ml of insecticide solution and one control bottle was coated with absolute ethanol.Alphacypermethrin and deltamethrin coated bottles were covered to keep them protected from light and to allow evaporation of the solvent overnight while the bottles coated with pirimiphos-methyl were left covered and protected from light for 8 hours.
Approximately 10 to 25 mosquitoes aged 4-5 days and fed with 10% honey solution were released in each bottle using an aspirator.The exposure time was set to 30 min, after which the mosquitoes were removed from the bottles and sorted into "alive" and "knocked-down" groups.Mosquitoes alive after 30 minutes of exposure to insecticide were considered resistant [47,48].Mosquitoes from the control bottle were labeled as unexposed.

Mosquito species identi cation
The legs of the resistant and unexposed mosquitoes were removed, and the rest of the body was immediately stored in RNA later.DNA extraction from the legs was performed following the protocol of Myriam and Cecile (2003) using the cetyltrimethyl ammonium bromide (CTAB) technique.Mosquito species identif ication was carried out to distinguish members of the An.gambiae s.l species complex using species-speci c PCR primers for An.gambiae s.s, and An.arabiensis.The PCR mix had a total volume of 20 µl composed of 5 µl of DNA, 8.4 µl of H 2 O, 2.5 µl of reaction buffer 1x, 0.1 µl of Taq DNA polymerase, 1 µl of MgCl 2 , and 1 µl of each primer (Universal forward primer: 5'-GTGTGCCCCTTCCTCGATGT-3'; An. gambiae s.s reverse primer:5'-CTGGTTTGGTCGGCACGTTT-3' and An.arabiensis reverse primer: 5'-AAGTGTCCTTCTCCATCCTA).The PCR cycling conditions were: 94°C for 5 min, followed by 25 ampli cation cycles (94°C for 30 s, 72°C for 30 s) and a nal elongation step at 72°C for 5 min.The PCR products were visualized using 1.5% agarose gels stained with 5 µl of BET.The DNA bands on the gel (390 bp for An.gambiae s.s. and 315 bp for An.arabiensis) were compared to a 1 kb reference ladder.
RNA extraction, RNA-Seq library preparation and sequencing Four-to-ve-day-old adult nonblood-fed female mosquitoes from the susceptible Kisumu strain were killed by freezing and stored at -80°C until RNA extraction.Mosquitoes from the eld populations unexposed to insecticides and those that survived to the highest dose of each insecticide were stored in RNAlater and shipped to the Entomology Branch laboratory at the CDC, Atlanta, USA, for RNA extraction, library preparation and sequencing.
RNA extraction was conducted using the Arcturus PicoPure RNA isolation kit (Life Technologies, USA) according to the manufacturer's instructions from three biological replicates with pools of 6 mosquitoes each from the following groups: mosquitoes from Djougou phenotyped as resistant to pirimiphosmethyl (DP), alphacypermethrin (DA) or deltamethrin (DD); mosquitoes from Djougou unexposed to insecticides (DU); mosquitoes from Bassila phenotyped as resistant to pirimiphos-methyl (BP) or deltamethrin (BD); mosquitoes from Bassila unexposed to insecticides (BU) and mosquitoes from the susceptible Kisumu strain (KIS).The Agilent 4200 TapeStation was used to measure the RNA concentration and integrity.Ribosomal RNA was depleted using the Ribo-Zero™ Magnetic Core Kit and Ribo-Zero™ rRNA Removal kit (Illumina, USA).Library preparation was carried out using the ScriptSeq v2 RNA-Seq Library Preparation Kit (Epicenter, Illumina) according to the manufacturer's instructions.Libraries were puri ed using Agencourt AMPure XP beads (Beckman Coulter, USA).The quantity and size distributions of the libraries were assessed using the Agilent DNA ScreenTape assay.Equimolar amounts of each library were pooled and sequenced (2x125 bp paired-end) on an Illumina HiSeq 2500 sequencer, using v2 chemistry.Sequencing was performed at the Biotechnology Core Facility at CDC, Atlanta, USA.

RNA-Seq data quality control mapping
Quality control was performed on the raw dmultiplexed paired end sequencing reads obtained from the sequencing center using FastqC [49].Reads from lane 1 and lane 2 were concatenated and trimmed to remove polyG tails, polyX at the 3' ends, bases that did not meet the minimum quality score of 20 and paired reads where one or both were shorter than 50 bp using Fastp [50].The trimmed reads (R1/R2) were aligned to the An.gambiae PEST reference genome (GenBank assembly identi er = GCA_000005575.2),directly downloaded from Vectorbase (release48) using 'subjunc', which is part of the subread aligner package, version 2.0.1 [51] with default settings.The alignments were ltered to remove low quality mapping reads (< 10) and sorted using SAMtools (v.1.10)[52].
FeatureCounts from the subread package version 2.0.1 [51,53] was used to count tags that overlapped the coding sequence (CDS) features by at least 1 bp in the sense orientation of the gene set AgambiaePest (structural annotation version = AgamP4.13).Tag counts represent the number of sequence reads that originated from a particular gene.The higher the number of counts, the more reads associated with that gene, and the assumption that there is a higher level of expression of that gene in the sample.The tag count result was used as input for differential expression analysis using edgeR, version 4.1.1[54] .
Differential expression analysis was performed using the package EdgeR to decide whether, for a given gene, an observed difference in tag counts was signi cantly different from what would be expected due to natural random variation.Thus, to remove the effect of noise and genes with very low expression, only genes where at least a tag count of 50 or more was obtained across all libraries were considered.
CalcNormFactors, a function of the edgeR package that uses the TMM (trimmed mean M-values) method, was used to normalize the number of tags between samples, by nding a set of scaling factors for the library sizes that minimize log-fold changes between samples for most genes.DEGs between the different comparisons were selected after multiple testing using the decideTests function, of the limma package [55].An absolute fold-change > 2 and FDR (false discovery Rate) ≤ 0.01 were used as statistical cutoffs to tag a gene as a DEG.
Comparisons were made between i) resistant eld mosquitoes and laboratory susceptible mosquitoes (R-S): (DA vs KIS; DD vs KIS; DP vs KIS; BD vs KIS; BP vs KIS); ii) unexposed eld mosquitoes and laboratory susceptible mosquitoes (C-S): (DU vs KIS; BU vs KIS) and iii) resistant eld mosquitoes and unexposed eld mosquitoes (R-C): (DA vs DU; DD vs DU; DP vs DU; BD vs BU; BP vs BU).R-S and C-S comparisons were made following the assumptions that constitutive resistance genes would be differentially expressed between both bioassay survivors and the unexposed mosquitoes when compared to the susceptible strain.The R-C comparison was made to detect gene expression that may be induced due to insecticide exposure.In addition, insecticidespeci c comparisons were made (DD vs BD and DP vs BP) to detect DEGs associated with the same insecticide across the different mosquito populations.
Detection of non-synonymous target site mutations site mutations detected from the variant calling analysis.Sorted bam les were used to count the read coverage at each genomic position and provide genotype likelihoods using the 'mpileup' methods of BCFtools [52,56] .Single nucleotide variants (SNVs) were then detected using bcftool 'call' and the resulting vcf le was ltered with vcfutils.pl[52,56] .

Gene ontology enrichment analysis
Gene ontology enrichment analysis was used to classify the differentially expressed genes based on the speci c biological functions using Goatools [57] .
This approach allows enrichment analysis to indicate which molecular functions, biological processes or cellular components were overrepresented (enriched) in a DEG list compared to an annotated list of the whole genome of An. gambiae obtained from blast2go [58] (Additional le 12) .The, GO term enrichment analysis of up-and downregulated genes was carried out using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) [59,60].The P-values used to assess signi cantly enriched GO terms associated with the DEG list were calculated based on Fisher's exact test and corrected by the Benjamini-Hochberg multiple testing correction method.Thus, an FDR adjusted P-value < 0.05 was used to report signi cantly enriched GO terms from the list of DEGs (Additional le 10).Visualization of the enriched GO terms was performed using the package 'ggplot2' following the protocol of Bonnot and others [61] .

RNAseq data validation using quantitative PCR
To validate the RNA-seq data, a set of six genes (SG7, CYP9K1, CYP6P3, COEJHE5E, RPS7, Actinc5) was used.RNA from three replicates of samples resistant to alphacypermethrin, deltamethrin or pirimiphos-methyl from Djougou (same batch of mosquitoes used for RNA sequencing) was used to synthesize cDNA using the HighCapacity cDNA reverse transcription kit (Applied Biosystems) with oligo-dT20 (NEB), according to the manufacturer's instructions.The primers used are listed in Additional le 13.Standard curves of Ct values for each gene were generated using a serial dilution of cDNA, allowing assessment of PCR e ciency.qPCR ampli cation was carried out on an Agilent Technologies Stratagene Mx3005P using PowerUp SYBR Green Master Mix (Applied Biosystems).cDNA from each sample was used as a template in a three-step program: 50°C for 2 minutes denaturation at 95°C for 10 minutes, followed by 40 cycles of 15 seconds at 95°C, 1 minute at 60°C and a nal step of 15 seconds at 95°C, 1 minute at 60°C, and 15 seconds at 95°C.The relative expression level and fold change (FC) of each target gene from resistant eld samples relative to the susceptible lab strain were calculated using the 2 − ΔΔCT method [62] incorporated in Python script (https://github.com/dany-gaga/from_Ct_to_logFoldChange).Housekeeping genes encoding ribosomal protein S17 (RPS17; AGAP010592) and Actin5C (AGAP000651) were used for normalization.The average mortalities of mosquitoes exposed to alphacypermethrin, deltamethrin and primiphos-methyl for 30 minutes are shown as percentages on the y-axis with 95% con dence intervals.Venn diagrams showing differentially expressed genes among resistant, susceptible and unexposed mosquito populations.
Panel A represents genes differentially expressed in Bassila and Djougou mosquitoes when exposed to deltamethrin; B represents genes differentially expressed when mosquitoes from both sites were exposed to pirimiphos-methyl; C represents genes differentially expressed when mosquitoes from both sites were exposed to pyrethroids.Each Venn diagram section shows the number of differentially expressed genes meeting each set of conditions and the P-values were adjusted for multiple testing based on FDR<0.01 and FC > 2.

1
AbbreviationsIR Insecticide resistance ITNs Insecticide-treated nets IRS indoor residual spraying DGE differential gene expression FC Fold-change CPs Cuticular proteins CYPs Cytochrome P450 genes GSTs Glutathione-S-transferases COE Carboxylesterases SNPs Single nucleotide polymorphisms GOE Gene ontology BP Biological process MF Molecular function CC Cellular component NMCP National Malaria Control Program LLINs Long-lasting insecticidal nets SGPs Salivary gland proteins RPS17 Ribosomal protein S1 CDS coding sequence GO Gene ontology

Figure 2 Gene
Figure 2

Figure 3 Key
Figure 3

Figure 7 Comparison
Figure 7

Table 1
Differential gene expression summary for alphacypermethrin, deltamethrin and pirimiphos-methyl

Table 2
Signi cantly differentially expressed genes of interest in R-S and C-S (FDR < 0.01 and FC > 2)